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Abstract. The augmented Lagrangiam method (ALM), widely used in quantum chemistry constrained 
optimization problems, is applied in the context of the nuclear Density Functional Theory (DFT) in the 
self-consistent constrained Skyrme Hartree-Fock-Bogoliubov (CHFB) variant. The ALM allows precise cal- 
culations of multidimensional energy surfaces in the space of collective coordinates that are needed to, e.g., 
determine fission pathways and saddle points; it improves accuracy of computed derivatives with respect 
to collective variables that are used to determine collective inertia; and is well adapted to supercomputer 
applications. 



PACS. 02.60.Pn Numerical optimization 
Functional Theory and extensions 



31.15.E Density-functional theory - 21.60.Jz Nuclear Density 



(N 
> 

m 



o 
o 



X 



1 Introduction 

The HFB equations of the superconducting nuclear DFT 
[T] can be viewed as a constrained nonlinear optimization 
problem in which the total energy of the nucleus, repre- 
sented by a functional of one-body densities, is minimized 
subject to constraints on the values of several independent 
variables. In addition to the usually imposed conditions on 
the number of particles (protons and neutrons), one is of- 
ten interested in constraining angular momentum compo- 
nents (to study nuclear behavior at nonzero angular mo- 
mentum) or nuclear multipole moments (or deformations) 
- to investigate the large amplitude collective motion, such 
as shape coexistence, fission or fusion. 

Constrained calculations are also used when going be- 
yond the standard single- reference DFT, e.g., within the 
Generator Coordinate Method [5] of the multi-reference 
DFT [2111] , where the constrained HFB solutions are used 
to generate a set of basis wave functions employed in fur- 
ther optimization. Another set of applications concerns 
the adiabatic approximation to the time-dependent HFB 
(ATDHFB) Qj^iOJjiSMl wherein derivatives with respect 
to collective coordinates are often approximated by finite- 
difference expressions [TU] . 

The fission problem is of particular interest as it in- 
volves many constrained calculations along collective de- 
grees of freedom representing families of mean fields char- 
acterizing fission pathways and nuclear dynamics during 
the fission process. In particular, care should be taken to 



identify saddle points in a multidimensional energy sur- 
face [TTlfT^fT^ . In this respect, constrained calculations 
in many variables can be very helpful as they can sepa- 
rate potential energy sheets that overlap when studied in 
reduced-deformation spaces [T4] . 

An effective approach to satisfy constraints is the method 
of Lagrange multipliers |15| . For example, the minimiza- 
tion of energy E at the condition that the nuclear quadrupole 
moment Q has an expectation value 5°, is equivalent to 
minimization of the Lagrangian function (or Routhian) 
E' = E + X{{Q) — g"), where the Lagrange multiplier A is 
determined from the condition 



(Q) = 9°. 



(1) 



In many cases, however, the procedure based on a 
linear constraint method (LCM) fails and the standard 
technique adopted is the method of quadratic constraint 
(the quadratic penalty method, QPM) fTBllTTlfTS] . In the 
above example, the corresponding Lagrangian function 
reads E' = E + c{{Q) — q^Y ■ As noted in early nuclear 
self-consistent applications P^[^ . results of calculations 
based on QPM strongly depend on the magnitude of c. 
For example, when the value of c is too small, one ends up 
with a solution having the constrained moment quite far 
away from the requested value. Increasing the value of c 
is often impossible as the self-consistent procedure ceases 
to converge. This is a serious deficiency of the method as 
it leaves important domains of the collective space unre- 
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solved thus obstructing (or even preventing) the theoreti- 
cal description. 

An effective procedure that avoids some of the difficul- 
ties pertaining to the standard LCM but not introducing 
the penalty term is the method proposed in [3T] and used 
in early CHF calculations of Refs. j^l^. and in CHFB 
calculations of Refs. [Mll^ . in which the A is changed it- 
eratively to satisfy the condition ([IJ at each step. Also, 
the constrained optimization problem is often treated by 
means of the conjugate gradient method [26j. However, ex- 
cept for a few cases [37] , most of the existing HFB solvers 
are based on a direct diagonalization approach and a mix- 
ing of intermediate solutions during the iteration process 
using a linear or Broyden mixing [28] . 

In this study, we demonstrate that the augmented La- 
grangian method {^[50] is an excellent alternative for 
nuclear-constrained HFB calculations. We show that the 
method always yields self-consistent solutions correspond- 
ing to requested values of constraints, independently of 
the value of the Lagrange multiplier selected. In this way, 
adopting the ALM, one can always access any region of the 
multi dimensional energy surface requested by the partic- 
ular physical phenomena investigated. At the same time, 
practical implementations of ALM do not require more 
computational resources as compared to QPM. A proce- 
dure, based on QPM but introducing a modifications of (j° 
during the iterations through a linear constraint has been 
used in [31]. While not based on the ALM algoritm, the 
spirit of this method is close to ALM. 

This paper is organized as follows. The method of La- 
grange multipliers, in its linear and quadratic variants, is 
briefly discussed in Sec. H] together with the augmented 
Lagrangian method. The ALM algorithm adopted in this 
work for diagonalization-based HFB solvers is laid out in 
Sec. |3l and the illustrative results are given in Sec. 01 Fi- 
nally, conclusions are contained in Sec. (5] 



2 The Method of Lagrange Multipliers 

A constrained optimization problem is usually specified in 
terms of equality and inequality constraints [TOllTTlfTB] . We 
consider here a finite-dimensional, equality-constrained non- 
linear optimization problem (ECP) of the form 

{mm£(x) 
(2) 
subject to gi{x) ^ qf, i = 1, 2, . . . , m, 

where we assume that £ : R" R (an objective func- 
tion) and gi : M" — > K (the constraint functions) are 
smooth functions, and n > m. The Lagrangian function 
E' : R"+" -> R associated with ECP is defined as 

rn 

E\x,\)^£{x) + J2^d9'^(x)-1°] 

i=l 

^£{x) + X'^[g{x)-q% (3) 
where A={Ai} is the vector of Lagrange multipliers. 



The following set of necessary and sufficient conditions 
allow the problem ([2]) to be formulated in terms of the 
Lagrangian function 

— The first-order necessary {local zero-slope) condition: 
if X* e R" is a local solution of ECP then there 
exists a unique vector A* S R™ such that (a;*, A*) is 
a stationary point for the Lagrangian function 

V,£;'(a;*,A*) = 0. (4) 

— The second-order necessary (local convexity) condition: 
if the functions £{x) and gi{x) are twice continuously 
differentiable, then the Hessian matrix of the Lagrangian 
function (with respect to x) must be positive semidef- 
inite at (a;*, A*) 

i?(s')(a;*,A*) >r 0. (5) 

In this context it should be mentioned, that Eq. ([5]) is 
the necessary and sufficient condition for convexity of 
the E (•, A*) function on its convex domain. 

— The second-order sufficient conditions: assume that £{x) 
and gi{x) are twice continuously differentiable, and let 
X* e R" and A* e R™ satisfy the equations 

V,e'{x*,\*) = 0, VxE\x*,\*)=0, (6) 

and the Hessian matrix H{E ){x, A) is positive definite 
at (cc*. A*) 

H{e'){x*,\*)>~0, (7) 
then the vector x* is a strict local solution of ECP 

In reality, these conditions are not always satisfied. For 
example, even when Eq. ([2]) has a solution, the Lagrangian 
function E could be unbounded [5^155] since the Hessian 
of the Lagrangian function is not necessarily positively 
defined. However, it can be shown that if ECP is con- 
vex^ i.e., £{x) is a convex (U-shaped) function and g{x) is 
an affine function, then the local constrained minimum is 
unique and represents the global minimum. 

Suppose that x{q) and A(q) are continuously differ- 
entiable functions giving the local minima of a family of 
ECPs (0) parameterized by q G R™, and x{q°) = x* , 
A(q°) = A*. This implies [H] that 

V,£[x{q)] = -\{q) or VE{q)^-\{q), (8) 

where the primal function E{q) is given by 

E{q) = £[xiq)] = min £ix). (9) 
g{x)=q 

The interpretation of the Lagrangian multipliers is there- 
fore that — A(<7) defines the parametrical gradient of E{q). 
When £{x) is the convex function, then a primal function 
E{q) turns out to be a convex function as well [Mj . 

The following sections illustrate the method of La- 
grangian multipliers for a one-dimensional problem. We 
first discuss the linear constraint method, then quadratic 
penalty method, and augmented Lagrangian method. 
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2.1 The linear constraint method 

The curve E{q) in Fig. [IJa) represents schematically how 
the primal function might vary as a function of q. It should 
be noted that we analyze the primal function E{q), not 
the objective function £{x). 




Fig. 1. (Color online) (a) Geometric interpretation of the lin- 
ear constraint method (LCM) . This method is not applicable to 
the shaded region, where the primal function E{q) is concave. 

(b) Behavior of penalty function for one-dimensional QPM. 

(c) Geometric interpretation of the ALM iterations (based on 
Ref. [IT])- See text for details. 



To solve the optimization problem, one can draw the 
tangent to the function E{q) at the point q = q^ , and use 
this line as the abscissa axis of a new coordinate system 
rotated by the angle a" with tana" = — A" = ^(9°), see 
Eq. . The ordinate axis in the new frame corresponds to 
the Lagrangian function E {q, A") = E{q^) — X°{q — q°). An 
unconstrained minimization of E gives the local minimum 
in the rotated system at the requested point q'^. In this 
way, the constrained minimization of E(q) is achieved by 
an unconstrained minimization of E {q, X^). 

However, as discussed above, the constrained mini- 
mization procedure with LCM can be applied only in 
the regions of E{q) which are convex. In the concave (fl- 
shaped) region, shaded in Fig.[TJa), the function E in the 
rotated frame has a maximum at point q^ . The minimiza- 
tion procedure docs not yield a stable solution around the 
maximum, and the constrained calculation with a linear 
constraint function fails to converge in the whole shaded 
region, see discussion in Refs. P^I^I^I] and in Sec. 7.6 of 
Ref. [2|. 

2.2 The quadratic constraint approach 

The local convexity assumption ([5]) plays a crucial part 
in solving the constrained problem The QPM can be 
applied when the convexity of original ECP is not pre- 
served. In such cases, one approximates the original con- 
strained problem by an unconstrained minimization prob- 
lem that involves a penalty for violation of the constraints, 
sec Refs. [16l[l7l[18] , and also Refs. [2j[19l[20]: 

min£:^(a;) = min \ £{x) + c V[5i(a;) - q°]^ \ 

I 1=1 ) 

= mm{£ix) + c\\g{x)-qY}, (10) 

where c > is called penalty parameter and || • || denotes 
Euclidean norm. It should be noted that if c is taken suf- 
ficiently large, then the local convexity condition can be 
shown to hold for the Lagrangian function E^{x). 

Figure [IJb) shows the same one-dimensional case as in 
Fig. (Ha), but for the QPM. The primal function E{q) is 
plotted with a solid (black) line, while the penalty function 
is plotted with dashed lines around the requested point q'^ 
for two different values ci and C2 of the penalty parameter 
c. The resulting Lagrangian functions E^{q) = E{q)+c{q— 
q^Y are indicated. 

It is immediately seen in Fig. [TJb) that the minimum 
of the Lagrangian function does not correspond to but 
rather to the values qi and q2 corresponding to the penalty 
parameter ci and C2 , respectively. One can obtain the val- 
ues of the function E{q) in a broad range by changing the 
requested point g° or the penalty parameter c, or both. 
But one can neither predict in advance which value q will 
be reached, nor to produce a regular mesh of values, which 
is often of interest. 

We thus see that the main drawback of the QPM is 
that it never delivers exactly the requested constraint val- 
ues. If one denotes the solution to unconstrained problem 
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(dni) by x*{c) (i.e., E'^{x*{c)) w £{x*) = E{q^)), then it 
has been shown \JT that lim x*{c) = x* . However, the 

c— ^oo 

Hessian matrix H{E^){x) is ill-defined for large values c. 
One is therefore forced to make a compromise between 
satisfying the constraints and having a well-conditioned 
problem when using QPM pni55] . 



One can see in Fig. [TJc) that if A*^ is sufficiently close to 
A* and/or c'^ is sufficiently large, the next multiplier A'^^^ 
will be closer to A* than A*^ (see also Ch. 4 of [17] and 
Ch. 17 of [IH|). The general iterative algorithm, including 
an adjustment of c^, can be found in Noccdal and Wright 



2.3 The augmented Lagrangian method 

The ALM can be viewed as a combination of LCM and 
QPM. It was introduced as a computational tool in 1969, 
by Hestenes [12] and Powell ^U\, as an attempt to solve 
the difficulties of linear and quadratic constraint methods. 

Let us begin by introducing the augmented Lagrangian 
function 

Elix, A) = £{x) + \^[g{x) - q"] + c\\g{x) - q"f , (11) 
which is the Lagrangian function for the problem: 

min{£ix)+c\\gix)-qY} 

^ (12) 
subject to g{x) ~ , 

that has the same local minima as original ECP ©. The 
gradient of E^{x, A) with respect to x is 

V,£;l(a;, A) = V£{x) + Vg{x){\ + 2c[g{x) - q°]} 

^V^e'{x,\), (13) 

where 

~X = X + 2c[g{x)-q"]. (14) 

If A*"' is a good approximation to the solution A*, then 
it is possible to approach the optimum x* through the 
unconstrained minimization of E^f.{-,\'') without using 
large values of the corresponding penalty constant c'" . The 
only condition is that c'^ is sufficiently large to ensure that 
the augmented Lagrangian function E^f, is locally convex 
with respect to x. 

It has been shown in Refs. [2^[5IJ| that the use of the 
iterative Lagrange multipliers 

^k+i ^x'' + 2c''[g{x'') - q°] (15) 

leads to x'^ which minimize E^^ ( •, A'^). 

Figure djc) provides a geometric interpretation of the 
iteration ([TS]) . To understand this figure, note that if x'^ 
minimizes £"^4. ( ■, A*"'), then the vector q'^ = g{x^) mini- 
mizes E{q) + [\''Y[q - q"] + c'=j|q - qOf . Hence 

V{E{q) + c'^\\q-qY]\^^^.= 

= Vi?;.(q^) = -A^ (16) 

and 

VE{q'') = - (A'^ + 2c'=[q'= - q°]) = 

= -(A* + 2c*[g(a;^)-q"]) =-A'=+i. (17) 



3 The ALM Algorithm 

This section provides guidance on how to implement the 
ALM given a working QPM algorithm (which is the stan- 
dard way of carrying out constrained minimization with 
HFB solvers based on diagonalization iterations) . 

Let us consider, for simplicity, a solver which defines 
the energy of the system E{q), and q is the expectation 
value of the (quadruple) operator Q, i.e., q = (Q). If one 
wishes to compute E{q^) for a nuclear state with a re- 
quested quadruple deformation q^ , the minimized quan- 
tity in QPM is 

E'{q) = E{q)+c{q-qy. (18) 

Taking the variational derivative of (|18|) . the resulting 
mean-field potential can be written as 

h' = h + 2c{q-q°)Q. (19) 

In our study the penalty parameter c we keep fixed. With 
an appropriate value of c, the self-consistent procedure 
yields a solution with quadrupole deformation, which is 
close to the requested value q^ . 

Let us now apply the ALM. By taking the variational 
derivative of 

E'{q) = E{q) + \{q - q°) + c{q - q^)\ (20) 

the resulting mean-field potential becomes 

h' = h + 2c{q-q'^{\))Q, (21) 

where q^{\) ^ q" ~ A/2c. 

Comparing Eqs. ([T^ and ([?T|) . we see that the use of 
ALM practically does not need any change in the part of 
the solver that deals with constrains. One simply needs to 
substitute q° g°(A). The new information that needs to 
be supplied is the way A is updated during the iteration 
process. According to Eq. ([15]), the value A'^'*'"'^ depends on 
the previous value A*^' as: 

A^+i = A^ + 2c(q-g"). (22) 

The iterations can start from a zero value. A" = 0. This 
is a well-defined starting point as for A = 0, ALM reduces 
itself to QPM for which one assumes the solver is already 
working. 

Comparing Eqs. (|19p and ([HJ, one can see that in 
ALM the originally requested value is replaced by an 
effective value 9" (A) which is adjusted during the iteration 
process according to Eq. (|22p . At the end of iterations, of 
course, one ends up with the originally requested value q^ . 
The generalization of the ALM algorithm to the case of 
many constraint variables is straightforward. 



A. Staszczak et al.: Augmented Lagrangian Method for Constrained Nuclear Density Functional Theory 



5 



4 Results 

The ALM algorithm of Sec. [3] has been implemented and 
tested with the HFB solvers HFODD v.2.43c ^ and HF- 
BTHO |2Z] • The illustrative examples of calculations pre- 
sented below concern the spontaneous fission of ^^^Fm. 
We used the solver HFODD, which is capable of treating 
simultaneously all the possible collective degrees of free- 
dom that might appear on the way to fission. 

In the particle-hole channel, we employed the SkM* 
energy density functional [38]. In the pairing channel, we 
adopted a seniority pairing force with the strength param- 
eters fitted to reproduce the experimental gaps in ^^^Fm 
[55] . The single-particle basis consisted of the lowest 1,140 
stretched states originating from the lowest 31 major oscil- 
lator shells. The details of the calculations can be found in 
Ref . [13] . We wish to remark only that special care should 
be taken when combining ALM with the Broyden mixing 
[28]. 

Figured] shows the results of constrained DFT calcula- 
tions on a two-dimensional Cartesian grid of quadrupole 
(Q20) and octupole {Q30) moments. The requested val- 
ues of constrained moments correspond to the grid points. 
While the ALM yields these values very precisely, the stan- 
dard QPM yields multipole moments that are often very 
different from the requested ones. In particular, QPM is 
unable to cover the whole area of interest. 



T ' 1 ' 1 ' 1 ' 1 ' r 




I 1 I 1 I 1 I 1 I 1 I 

100 120 140 160 180 200 

Quadrupole moment Q (b) 



Fig. 2. (Color online) A comparison between the ALM (black 
squares) and the standard QPM (open squares) for the con- 
strained self-consistent convergence scheme. The HFB calcu- 
lations were carried out for the total energy surface of ^®^Fm 
in a two-dimensional plane of elongation, Q20, and reflection- 
asymmetry, Qao . Although QPM often fails to produce a solu- 
tion at the required values of constrained variables (Q20, Q30) 
on a rectangular grid, ALM performs very well in all cases. 

Figure [3] shows the total energy surface of ^^^Fm in 
the Q20-Q30 plane obtained with ALM. Comparing with 
Fig. 121 one can see interesting physics in the region which 
was inaccessible by QPM, namely the appearance of the 



second (fusion) valley at large values of Q30 separated 
from the spontaneous fission valley by a steep ridge. 




Quadrupole moment Q20 (b) 



Fig. 3. (Color online) Two-dimensional total energy surface for 
^^■^Fm calculated with SkM* energy density functional using 
the ALM in the plane of collective coordinates Q20-Q30- The 
fission and fusion pathways are marked. The difference between 
contour lines is 5 MeV. 



5 Conclusions 

The augmented Lagrangiam method to solve a constrained 
nonlinear problem of CHFB has been compared with the 
standard variants of the method of Lagrange multipli- 
ers, i.e., quadratic penalty method and linear constraint 
method. We discuss the numerical strategy beyond QPM 
and ALM algorithms and show how to implement ALM 
in HFB solvers based on the diagonalization approach. 

Compared to QPM, we find ALM to be superior: it 
enables precise constrained calculations in many dimen- 
sions thus uncovering regions of collective space that are 
not accessible with the standard method. The method is 
well adapted to supercomputer applications and its stabil- 
ity makes it a tool of choice for large-scale CHFB calcu- 
lations, such as computations of multidimensional fission 
pathways discussed in this work. 
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